home *** CD-ROM | disk | FTP | other *** search
/ Enter 2004 January / enter-2004-01.iso / files / maxima-5.9.0.exe / {app} / share / maxima / 5.9.0 / src / numerical / slatec / zrati.lisp < prev    next >
Encoding:
Text File  |  2003-02-09  |  5.6 KB  |  156 lines
  1. ;;; Compiled by f2cl version 2.0 beta 2002-05-06
  2. ;;; 
  3. ;;; Options: ((:prune-labels nil) (:auto-save t) (:relaxed-array-decls t)
  4. ;;;           (:coerce-assigns :as-needed) (:array-type ':simple-array)
  5. ;;;           (:array-slicing nil) (:declare-common nil)
  6. ;;;           (:float-format double-float))
  7.  
  8. (in-package "SLATEC")
  9.  
  10.  
  11. (let ((czeror 0.0)
  12.       (czeroi 0.0)
  13.       (coner 1.0)
  14.       (conei 0.0)
  15.       (rt2 1.4142135623730951))
  16.   (declare (type double-float rt2 conei coner czeroi czeror))
  17.   (defun zrati (zr zi fnu n cyr cyi tol)
  18.     (declare (type f2cl-lib:integer4 n)
  19.              (type (simple-array double-float (*)) cyr cyi)
  20.              (type double-float zr zi fnu tol))
  21.     (prog ((i 0) (id 0) (idnu 0) (inu 0) (itime 0) (k 0) (kk 0) (magz 0)
  22.            (ak 0.0) (amagz 0.0) (ap1 0.0) (ap2 0.0) (arg 0.0) (az 0.0)
  23.            (cdfnui 0.0) (cdfnur 0.0) (dfnu 0.0) (fdnu 0.0) (flam 0.0)
  24.            (fnup 0.0) (pti 0.0) (ptr 0.0) (p1i 0.0) (p1r 0.0) (p2i 0.0)
  25.            (p2r 0.0) (rak 0.0) (rap1 0.0) (rho 0.0) (rzi 0.0) (rzr 0.0)
  26.            (test 0.0) (test1 0.0) (tti 0.0) (ttr 0.0) (t1i 0.0) (t1r 0.0))
  27.       (declare
  28.        (type double-float t1r t1i ttr tti test1 test rzr rzi rho rap1 rak p2r
  29.         p2i p1r p1i ptr pti fnup flam fdnu dfnu cdfnur cdfnui az arg ap2 ap1
  30.         amagz ak)
  31.        (type f2cl-lib:integer4 magz kk k itime inu idnu id i))
  32.       (setf az (zabs zr zi))
  33.       (setf inu (f2cl-lib:int fnu))
  34.       (setf idnu (f2cl-lib:int-sub (f2cl-lib:int-add inu n) 1))
  35.       (setf magz (f2cl-lib:int az))
  36.       (setf amagz
  37.               (coerce (the f2cl-lib:integer4 (f2cl-lib:int-add magz 1))
  38.                       'double-float))
  39.       (setf fdnu (coerce (the f2cl-lib:integer4 idnu) 'double-float))
  40.       (setf fnup (max amagz fdnu))
  41.       (setf id (f2cl-lib:int-sub idnu magz 1))
  42.       (setf itime 1)
  43.       (setf k 1)
  44.       (setf ptr (/ 1.0 az))
  45.       (setf rzr (* ptr (+ zr zr) ptr))
  46.       (setf rzi (* (- ptr) (+ zi zi) ptr))
  47.       (setf t1r (* rzr fnup))
  48.       (setf t1i (* rzi fnup))
  49.       (setf p2r (- t1r))
  50.       (setf p2i (- t1i))
  51.       (setf p1r coner)
  52.       (setf p1i conei)
  53.       (setf t1r (+ t1r rzr))
  54.       (setf t1i (+ t1i rzi))
  55.       (if (> id 0) (setf id 0))
  56.       (setf ap2 (zabs p2r p2i))
  57.       (setf ap1 (zabs p1r p1i))
  58.       (setf arg (/ (+ ap2 ap2) (* ap1 tol)))
  59.       (setf test1 (f2cl-lib:fsqrt arg))
  60.       (setf test test1)
  61.       (setf rap1 (/ 1.0 ap1))
  62.       (setf p1r (* p1r rap1))
  63.       (setf p1i (* p1i rap1))
  64.       (setf p2r (* p2r rap1))
  65.       (setf p2i (* p2i rap1))
  66.       (setf ap2 (* ap2 rap1))
  67.      label10
  68.       (setf k (f2cl-lib:int-add k 1))
  69.       (setf ap1 ap2)
  70.       (setf ptr p2r)
  71.       (setf pti p2i)
  72.       (setf p2r (- p1r (- (* t1r ptr) (* t1i pti))))
  73.       (setf p2i (- p1i (+ (* t1r pti) (* t1i ptr))))
  74.       (setf p1r ptr)
  75.       (setf p1i pti)
  76.       (setf t1r (+ t1r rzr))
  77.       (setf t1i (+ t1i rzi))
  78.       (setf ap2 (zabs p2r p2i))
  79.       (if (<= ap1 test) (go label10))
  80.       (if (= itime 2) (go label20))
  81.       (setf ak (* (zabs t1r t1i) 0.5))
  82.       (setf flam (+ ak (f2cl-lib:fsqrt (- (* ak ak) 1.0))))
  83.       (setf rho (min (/ ap2 ap1) flam))
  84.       (setf test (* test1 (f2cl-lib:fsqrt (/ rho (- (* rho rho) 1.0)))))
  85.       (setf itime 2)
  86.       (go label10)
  87.      label20
  88.       (setf kk (f2cl-lib:int-sub (f2cl-lib:int-add k 1) id))
  89.       (setf ak (coerce (the f2cl-lib:integer4 kk) 'double-float))
  90.       (setf t1r ak)
  91.       (setf t1i czeroi)
  92.       (setf dfnu (+ fnu (f2cl-lib:int-sub n 1)))
  93.       (setf p1r (/ 1.0 ap2))
  94.       (setf p1i czeroi)
  95.       (setf p2r czeror)
  96.       (setf p2i czeroi)
  97.       (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
  98.                     ((> i kk) nil)
  99.         (tagbody
  100.           (setf ptr p1r)
  101.           (setf pti p1i)
  102.           (setf rap1 (+ dfnu t1r))
  103.           (setf ttr (* rzr rap1))
  104.           (setf tti (* rzi rap1))
  105.           (setf p1r (+ (- (* ptr ttr) (* pti tti)) p2r))
  106.           (setf p1i (+ (* ptr tti) (* pti ttr) p2i))
  107.           (setf p2r ptr)
  108.           (setf p2i pti)
  109.           (setf t1r (- t1r coner))
  110.          label30))
  111.       (if (or (/= p1r czeror) (/= p1i czeroi)) (go label40))
  112.       (setf p1r tol)
  113.       (setf p1i tol)
  114.      label40
  115.       (multiple-value-bind
  116.           (var-0 var-1 var-2 var-3 var-4 var-5)
  117.           (zdiv p2r p2i p1r p1i (f2cl-lib:fref cyr (n) ((1 n)))
  118.            (f2cl-lib:fref cyi (n) ((1 n))))
  119.         (declare (ignore var-0 var-1 var-2 var-3))
  120.         (f2cl-lib:fset (f2cl-lib:fref cyr (n) ((1 n))) var-4)
  121.         (f2cl-lib:fset (f2cl-lib:fref cyi (n) ((1 n))) var-5))
  122.       (if (= n 1) (go end_label))
  123.       (setf k (f2cl-lib:int-sub n 1))
  124.       (setf ak (coerce (the f2cl-lib:integer4 k) 'double-float))
  125.       (setf t1r ak)
  126.       (setf t1i czeroi)
  127.       (setf cdfnur (* fnu rzr))
  128.       (setf cdfnui (* fnu rzi))
  129.       (f2cl-lib:fdo (i 2 (f2cl-lib:int-add i 1))
  130.                     ((> i n) nil)
  131.         (tagbody
  132.           (setf ptr
  133.                   (+ cdfnur
  134.                      (- (* t1r rzr) (* t1i rzi))
  135.                      (f2cl-lib:fref cyr ((f2cl-lib:int-add k 1)) ((1 n)))))
  136.           (setf pti
  137.                   (+ cdfnui
  138.                      (+ (* t1r rzi) (* t1i rzr))
  139.                      (f2cl-lib:fref cyi ((f2cl-lib:int-add k 1)) ((1 n)))))
  140.           (setf ak (zabs ptr pti))
  141.           (if (/= ak czeror) (go label50))
  142.           (setf ptr tol)
  143.           (setf pti tol)
  144.           (setf ak (* tol rt2))
  145.          label50
  146.           (setf rak (/ coner ak))
  147.           (f2cl-lib:fset (f2cl-lib:fref cyr (k) ((1 n))) (* rak ptr rak))
  148.           (f2cl-lib:fset (f2cl-lib:fref cyi (k) ((1 n))) (* (- rak) pti rak))
  149.           (setf t1r (- t1r coner))
  150.           (setf k (f2cl-lib:int-sub k 1))
  151.          label60))
  152.       (go end_label)
  153.      end_label
  154.       (return (values nil nil nil nil nil nil nil)))))
  155.  
  156.